package scu.maqiang.homogenization;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2Q41;
import scu.maqiang.fes.FES2Q81;
import scu.maqiang.mesh.Mesh2Q4;
import scu.maqiang.mesh.Mesh2Q8;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class PossionSOTSProblem2DQ {

	public static void main(String[] args) {
		long begin = System.nanoTime();
	    int N = 16;
	    double eps= 1.0 / N;
		Mesh2Q4 cellMesh = new Mesh2Q4().square2D(32, 32);
	    ScalarFunc regionCellFunc = (xy, label, param) -> {
            return (Math.abs(xy[0] - 0.5) < 0.25 && Math.abs(xy[1] - 0.5) < 0.25)? 1.0:-1.0;
        };
	    cellMesh.setDomainLabel(regionCellFunc, null, 10);
//	    cellMesh.displayElements();
	    FES2Q41 fsCell = new FES2Q41(cellMesh);
	    double[][] FOCS = new double[2][fsCell.GetNdof()];
		double[][] SOCS = new double[4][fsCell.GetNdof()];
		double[][] D = new double[2][2];
		SOTS2DQ481 ss = new SOTS2DQ481(fsCell, null);
		ScalarFunc cf = (xyz, label, param) -> label == 10? 0.001:1;
		
		ScalarFunc[] gcf = {cf};
		ss.KxxCellSolution(gcf, null, FOCS, D, SOCS);
	    
		System.out.println("min: " + MVO.min(FOCS[0]) + "\tmax: " + MVO.max(FOCS[0]));
		System.out.println("min: " + MVO.min(FOCS[1]) + "\tmax: " + MVO.max(FOCS[1]));
		System.out.println("Homogenized Kxx: ");
		System.out.println(MVO.toString(D));
		System.out.println("min: " + MVO.min(SOCS[0]) + "\tmax: " + MVO.max(SOCS[0]));
		System.out.println("min: " + MVO.min(SOCS[1]) + "\tmax: " + MVO.max(SOCS[1]));
		System.out.println("min: " + MVO.min(SOCS[2]) + "\tmax: " + MVO.max(SOCS[2]));
		System.out.println("min: " + MVO.min(SOCS[3]) + "\tmax: " + MVO.max(SOCS[3]));
		cellMesh.toTecplot("CellFOCSQ.dat", FOCS);
		cellMesh.toTecplot("CellSOCSQ.dat", SOCS);

		Mesh2Q8 homoMesh = new Mesh2Q8().square(64, 64);
		FES2Q81 fsHomo = new FES2Q81(homoMesh);
		ScalarFunc homoCoef = (xy, label, param) -> D[0][0];
		SRMatrix homoA = new SRMatrix(fsHomo.GetNdof());
		ScalarFunc[] MyHomoCoef = {homoCoef};
		fsHomo.assembleStiff(MyHomoCoef, null, BVPType.COMMON, homoA);
		double[] homoRHS = new double[fsHomo.GetNdof()];
		fsHomo.assembleSource(new double[] {10.0}, BVPType.COMMON, homoRHS);
		fsHomo.applyBC_MBN(homoA, Direct.X, 1, 2, 3, 4);
		fsHomo.applyBC_MBN(homoRHS, Direct.X, 0.0, 1, 2, 3, 4);
		NewIterSSolver homoSolver = new NewIterSSolver(homoA);
		double[] homoSol = new double[fsHomo.GetNdof()];
		homoSolver.PCGSSOR(homoRHS, homoSol, 1.5, 1);
		
		Mesh2Q4 fineMesh = new Mesh2Q4().square2D(256, 256);
		ScalarFunc regionFineFunc = (xy, label, param) -> {
            double xCell = xy[0] * N - (int)(xy[0] * N);
            double yCell = xy[1] * N - (int)(xy[1] * N);
            return (Math.abs(xCell - 0.5) < 0.25 && Math.abs(yCell - 0.5) < 0.25)? 1.0:-1.0;
        };
        fineMesh.setDomainLabel(regionFineFunc, null, 10);
        FES2Q41 fsFine = new FES2Q41(fineMesh);
        SRMatrix fineA = new SRMatrix(fsFine.GetNdof());
        fsFine.assembleStiff(new ScalarFunc[] {cf}, null, BVPType.COMMON, fineA);
        double[] fineRHS = new double[fsFine.GetNdof()];
        fsFine.assembleSource(new double[] {10.0}, BVPType.COMMON, fineRHS);
        fsFine.applyBC_MBN(fineA, Direct.X, 1, 2, 3, 4);
        fsFine.applyBC_MBN(fineRHS, Direct.X, 0.0, 1, 2, 3, 4);
        double[] fineSol = new double[fsFine.GetNdof()];
        NewIterSSolver fineSolver = new NewIterSSolver(fineA);
        fineSolver.PCGSSOR(fineRHS, fineSol, 1.5, 1);
        
        double[] HomoFineSol = new double[fsFine.GetNdof()];
        double[] FirstSol = new double[fsFine.GetNdof()];
        double[] SecondSol = new double[fsFine.GetNdof()];
        ss.assembleHeatSOTSSolution_Square(fsFine, HomoFineSol, FirstSol, SecondSol, fsHomo, homoSol, fsCell, FOCS, SOCS, N, eps);
        
        System.out.println("Homo Solution");
        System.out.println("Homo min: " + MVO.min(HomoFineSol) + "\tmax: " + MVO.max(HomoFineSol));
        System.out.println("First Order Approximation: ");
        System.out.println("FOA min: " + MVO.min(FirstSol) + "\tmax: " + MVO.max(FirstSol));
        System.out.println("Second Order Approximation: ");
        System.out.println("SOA min: " + MVO.min(SecondSol) + "\tmax: " + MVO.max(SecondSol));
        
        double[][] wholeSol = new double[][]{HomoFineSol, FirstSol, SecondSol, fineSol};
        fineMesh.toTecplot("SOTSSolutions.dat", wholeSol);
        double L2Norm = fsFine.computeL2Norm(fineSol);
        double[] error0 = MVO.add(fineSol, -1.0, HomoFineSol);
		double[] error1 = MVO.add(fineSol, -1.0, FirstSol);
		double[] error2 = MVO.add(fineSol, -1.0, SecondSol);
		System.out.println("L2RError:");
		System.out.println("Error 0: " + fsFine.computeL2Norm(error0) / L2Norm);
		System.out.println("Error 1: " + fsFine.computeL2Norm(error1) / L2Norm);
		System.out.println("Error 2: " + fsFine.computeL2Norm(error2) / L2Norm);
//		
		double H1SemiNorm = fsFine.computeH1SemiNorm(fineSol);
		System.out.println("H1SemiRError:");
		System.out.println("Error 0: " + fsFine.computeH1SemiNorm(error0) / H1SemiNorm);
		System.out.println("Error 1: " + fsFine.computeH1SemiNorm(error1) / H1SemiNorm);
		System.out.println("Error 2: " + fsFine.computeH1SemiNorm(error2) / H1SemiNorm);		
		
		System.out.println("LInftyRError");
		System.out.println("Error 0: " + MVO.maxA(error0) / MVO.maxA(fineSol));
		System.out.println("Error 1: " + MVO.maxA(error1) / MVO.maxA(fineSol));
		System.out.println("Error 2: " + MVO.maxA(error2) / MVO.maxA(fineSol));
		long end = System.nanoTime();
		System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
	}

}
